clear all; close all; clc;
sigma = 13;
N = 1e6;
y_t = 0:0.1:65;
p_t = y_t/sigma^2 .* ...
exp(- y_t.^2 / 2 / sigma^2 );
a = min(y_t); b = max(y_t);
maxp = max(p_t);
x1 = rand(1, N); x2 = rand(1, N);
x1r = x1*(b-a) + a; x2r = x2*maxp;
ind = find(x2r < x1r/sigma^2 .* ...
exp(- x1r.^2 / 2 / sigma^2));
y_n = x1r(ind);
[h, y] = hist(y_n, 30);
inte = sum(h) * (y(2) - y(1)); p = h / inte;
figure(1); plot(y_t, p_t, 'g', 'LineWidth', 5);
hold on; bar(y, p); hold off
xlabel('y'); ylabel('p(y)');
